﻿Namespace Imaging



    Public Class ColorSpaceConverter



#Region "Color Temperature"


        Protected Shared ReadOnly ReciprocalTemperature() As Double = {Double.MinValue, 0.00001, 0.00002, 0.00003, 0.00004, 0.00005, _
                                                                       0.00006, 0.00007, 0.00008, 0.00009, 0.0001, 0.000125, _
                                                                       0.00015, 0.000175, 0.0002, 0.000225, 0.00025, 0.000275, _
                                                                       0.0003, 0.000325, 0.00035, 0.000375, 0.0004, 0.000425, _
                                                                       0.00045, 0.000475, 0.0005, 0.000525, 0.00055, 0.000575, _
                                                                       0.0006}

        Protected Shared ReadOnly UVT(,) As Double = {{0.18006, 0.26352, -0.24341}, _
                                                      {0.18066, 0.26589, -0.25479}, _
                                                      {0.18133, 0.26846, -0.26876}, _
                                                      {0.18208, 0.27119, -0.28539}, _
                                                      {0.18293, 0.27407, -0.3047}, _
                                                      {0.18388, 0.27709, -0.32675}, _
                                                      {0.18494, 0.28021, -0.35156}, _
                                                      {0.18611, 0.28342, -0.37915}, _
                                                      {0.1874, 0.28668, -0.40955}, _
                                                      {0.1888, 0.28997, -0.44278}, _
                                                      {0.19032, 0.29326, -0.47888}, _
                                                      {0.19462, 0.30141, -0.58204}, _
                                                      {0.19962, 0.30921, -0.70471}, _
                                                      {0.20525, 0.31647, -0.84901}, _
                                                      {0.21142, 0.32312, -1.0182}, _
                                                      {0.21807, 0.32909, -1.2168}, _
                                                      {0.22511, 0.33439, -1.4512}, _
                                                      {0.23247, 0.33904, -1.7298}, _
                                                      {0.2401, 0.34308, -2.0637}, _
                                                      {0.24792, 0.34655, -2.4681}, _
                                                      {0.25591, 0.34951, -2.9641}, _
                                                      {0.264, 0.352, -3.5814}, _
                                                      {0.27218, 0.35407, -4.3633}, _
                                                      {0.28039, 0.35577, -5.3762}, _
                                                      {0.28863, 0.35714, -6.7262}, _
                                                      {0.29685, 0.35823, -8.5955}, _
                                                      {0.30505, 0.35907, -11.324}, _
                                                      {0.3132, 0.35968, -15.628}, _
                                                      {0.32129, 0.36011, -23.325}, _
                                                      {0.32931, 0.36038, -40.77}, _
                                                      {0.33724, 0.36051, -116.45}}

        ''' <summary>
        ''' Converts a color from XYZ color space to Correlated Color Temperature
        ''' </summary>
        ''' <param name="XYZ">Adjusted XYZ coordinates (0..1) of XYZ color sapce</param>
        ''' <returns>The corresponding color temperature between 1667 and 100000</returns>
        ''' <remarks>Returns 0 if input coordinates are out of range.</remarks>
        Public Shared Function XYZtoCorColorTemp(ByVal XYZ As XYZTriplet) As Integer
            Dim vU As Double
            Dim vV As Double
            Dim p As Double
            Dim di As Double
            Dim dm As Double
            Dim Index As Integer
            Dim result As Integer = -1

            ' protect against possible divide-by-zero failure
            If (XYZ.X >= 1.0E-20) AndAlso (XYZ.Y >= 1.0E-20) AndAlso (XYZ.Z >= 1.0E-20) Then
                vU = (4.0 * XYZ.X) / (XYZ.X + 15.0 * XYZ.Y + 3.0 * XYZ.Z)
                vV = (6.0 * XYZ.Y) / (XYZ.X + 15.0 * XYZ.Y + 3.0 * XYZ.Z)
                dm = 0.0

                For i As Integer = 0 To 30
                    di = (vV - UVT(i, 1)) - UVT(i, 2) * (vU - UVT(i, 0))
                    If ((i > 0) AndAlso (((di < 0.0) AndAlso (dm >= 0.0)) OrElse ((di >= 0.0) AndAlso (dm < 0.0)))) Then
                        Index = i
                        Exit For ' found lines bounding (us, vs) : i-1 and i
                    End If

                    dm = di
                Next
                If Index = -1 Then Return 0 ' bad XYZ input, color temp would be less than minimum of 1666.7 degrees, or too far towards blue

                di = di / Math.Sqrt(1.0 + UVT(Index, 2) * UVT(Index, 2))
                dm = dm / Math.Sqrt(1.0 + UVT(Index - 1, 2) * UVT(Index - 1, 2))
                p = dm / (dm - di)     ' p = interpolation parameter, 0.0 : i-1, 1.0 : i 
                result = Math.Round(1.0 / ((ReciprocalTemperature(Index) - ReciprocalTemperature(Index - 1)) * (p + ReciprocalTemperature(Index - 1))))
            End If
            Return result

        End Function
#End Region

        ''' <summary>Approximates the hue of a specific wavelength</summary>
        ''' <param name="Wavelength">A wavelength between 350nm and 780nm</param>
        ''' <returns>The approximated hue of the wavelength</returns>
        Public Shared Function Wavelength2Hue(ByVal Wavelength As Integer) As Double
            Dim Gamma As Double = 1.0
            Dim IntensityMax As Integer = 255
            Dim Blue As Double = 0.0
            Dim Green As Double = 0.0
            Dim Red As Double = 0.0
            Dim Factor As Double


            If (Wavelength >= 350 AndAlso Wavelength <= 439) Then
                Red = -(Wavelength - 440.0) / (440.0 - 350.0)
                Blue = 1.0
            ElseIf (Wavelength >= 440 AndAlso Wavelength <= 489) Then
                Green = (Wavelength - 440.0) / (490.0 - 440.0)
                Blue = 1.0
            ElseIf (Wavelength >= 490 AndAlso Wavelength <= 509) Then
                Green = 1.0
                Blue = -(Wavelength - 510.0) / (510.0 - 490.0)
            ElseIf (Wavelength >= 510 AndAlso Wavelength <= 579) Then
                Red = (Wavelength - 510.0) / (580.0 - 510.0)
                Green = 1.0
            ElseIf (Wavelength >= 580 AndAlso Wavelength <= 644) Then
                Red = 1.0
                Green = -(Wavelength - 645.0) / (645.0 - 580.0)
            ElseIf (Wavelength >= 645 AndAlso Wavelength <= 780) Then
                Red = 1.0
            End If


            If (Wavelength >= 350 AndAlso Wavelength <= 419) Then
                Factor = 0.3 + 0.7 * (Wavelength - 350.0) / (420.0 - 350.0)
            ElseIf (Wavelength >= 420 AndAlso Wavelength <= 700) Then
                Factor = 1.0
            ElseIf (Wavelength >= 701 AndAlso Wavelength <= 780) Then
                Factor = 0.3 + 0.7 * (780.0 - Wavelength) / (780.0 - 700.0)
            Else
                Factor = 0.0
            End If

            Dim R As Integer
            Dim G As Integer
            Dim B As Integer

            If Red = 0 Then R = 0 Else R = Math.Round(IntensityMax * Math.Pow(Red * Factor, Gamma))
            If Green = 0 Then G = 0 Else G = Math.Round(IntensityMax * Math.Pow(Green * Factor, Gamma))
            If Blue = 0 Then B = 0 Else B = Math.Round(IntensityMax * Math.Pow(Blue * Factor, Gamma))

            Return RGB2HSV(New RGBTriplet(Drawing.Color.FromArgb(R, G, B))).Hue

        End Function

        Public Shared Function Hue2Wavelength(ByVal Hue As Double) As Integer

        End Function

#Region "Color Space Conversion"

#Region "From RGB"

        ''' <summary>Converts floating point RGB to Radiance RGBE</summary>
        ''' <param name="RGB">Input RGB (values range from 0..1)</param>
        ''' <returns>RGBE (values in range of 0..255)</returns>
        Public Shared Function RGB2RGBE(ByVal RGB As RGBTriplet) As RGBEQuadlet
            Dim val As Double = Math.Min(Math.Max(Math.Max(RGB.Red, RGB.Green), RGB.Blue), 1)
            Dim exp As Integer
            Dim result As RGBEQuadlet

            If val >= 1.0E-32 Then
                val = ANSI_C_Math.frexp(val, exp) * 256.0 / val
                result.Red = Math.Min(Math.Max(RGB.Red * val, 0), 255)
                result.Green = Math.Min(Math.Max(RGB.Green * val, 0), 255)
                result.Blue = Math.Min(Math.Max(RGB.Blue * val, 0), 255)
                result.Exponent = exp + 128
            End If

            Return result
        End Function

        ''' <summary>Converts floating point RGB to Lab color space</summary>
        ''' <param name="RGB">Input RGB (values range from 0..1)</param>
        Public Shared Function RGB2Lab(ByVal RGB As RGBTriplet) As LabTriplet
            Return XYZ2Lab(RGB2XYZ(RGB))
        End Function

        ''' <summary>Converts floating point RGB to <see cref="Drawing.Color"/></summary>
        ''' <param name="RGB">Input RGB (values range from 0..1)</param>
        Public Shared Function RGB2Color(ByVal RGB As RGBTriplet, Optional ByVal Gamma As Double = 1) As Drawing.Color
            Return Drawing.Color.FromArgb(Math.Max(Math.Min(255, RGB.Red ^ (1 / Gamma) * 255), 0), _
                                          Math.Max(Math.Min(255, RGB.Green ^ (1 / Gamma) * 255), 0), _
                                          Math.Max(Math.Min(255, RGB.Blue ^ (1 / Gamma) * 255), 0))
        End Function

        ''' <summary>Converts floating point RGB to HSV color space</summary>
        ''' <param name="RGB">Input RGB (values range from 0..1)</param>
        Public Shared Function RGB2HSV(ByVal RGB As RGBTriplet) As HSVTriplet
            Dim vMax As Double = Math.Max(RGB.Red, Math.Max(RGB.Green, RGB.Blue))
            Dim dMinMax As Double = vMax - Math.Min(RGB.Red, Math.Min(RGB.Green, RGB.Blue))

            Dim HSV As HSVTriplet

            HSV.Value = vMax

            If (dMinMax = 0) Then
                HSV.Hue = 0
                HSV.Saturation = 0
            Else
                HSV.Saturation = dMinMax / vMax

                Dim dR As Double = (((vMax - RGB.Red) / 6) + (dMinMax / 2)) / dMinMax
                Dim dG As Double = (((vMax - RGB.Green) / 6) + (dMinMax / 2)) / dMinMax
                Dim dB As Double = (((vMax - RGB.Blue) / 6) + (dMinMax / 2)) / dMinMax

                If (RGB.Red = vMax) Then : HSV.Hue = dB - dG
                ElseIf (RGB.Green = vMax) Then : HSV.Hue = (1 / 3) + dR - dB
                ElseIf (RGB.Blue = vMax) Then : HSV.Hue = (2 / 3) + dG - dR
                End If

                If (HSV.Hue < 0) Then HSV.Hue += 1
                If (HSV.Hue > 1) Then HSV.Hue -= 1
            End If

            Return HSV
        End Function

        'TODO: consider WhitePoint/RGB Working Space
        'TODO: Gamma Correction

        ''' <summary>Converts floating point RGB to XYZ color space</summary>
        ''' <param name="RGB">Input RGB (values range from 0..1)</param>
        Public Shared Function RGB2XYZ(ByVal RGB As RGBTriplet) As XYZTriplet
            Dim vR As Double = RGB.Red
            Dim vG As Double = RGB.Green
            Dim vB As Double = RGB.Blue

            If (vR > 0.04045) Then vR = ((vR + 0.055) / 1.055) ^ 2.4 Else vR = vR / 12.92
            If (vG > 0.04045) Then vG = ((vG + 0.055) / 1.055) ^ 2.4 Else vG = vG / 12.92
            If (vB > 0.04045) Then vB = ((vB + 0.055) / 1.055) ^ 2.4 Else vB = vB / 12.92

            ' TODO: remove the following lines. Only valid for sRGB!
            Dim Result As XYZTriplet
            Result.X = vR * 0.4124564 + vG * 0.3575761 + vB * 0.1804375
            Result.Y = vR * 0.2126729 + vG * 0.7151522 + vB * 0.072175
            Result.Z = vR * 0.0193339 + vG * 0.119192 + vB * 0.9503041
            Return Result

        End Function

#End Region

#Region "From XYZ"


        Public Shared Function XYZ2RGB(ByVal XYZ As XYZTriplet) As RGBTriplet
            Return XYZ2RGB(XYZ.X, XYZ.Y, XYZ.Z)
        End Function

        'TODO: consider WhitePoint/RGB Working Space
        'see http://www.brucelindbloom.com/Eqn_RGB_XYZ_Matrix.html
        'and http://www.brucelindbloom.com/index.html?WorkingSpaceInfo.html#Specifications
        Public Shared Function XYZ2RGB(ByVal X As Double, ByVal Y As Double, ByVal Z As Double) As RGBTriplet

            Dim vX = Math.Max(Math.Min(X, 1), 0)
            Dim vY = Math.Max(Math.Min(Y, 1), 0)
            Dim vZ = Math.Max(Math.Min(Z, 1), 0)

            Dim vR As Double = vX * 3.2406 + vY * -1.5372 + vZ * -0.4986
            Dim vG As Double = vX * -0.9689 + vY * 1.8758 + vZ * 0.0415
            Dim vB As Double = vX * 0.0557 + vY * -0.204 + vZ * 1.057

            Dim Result As RGBTriplet

            If vR > 0.0031308 Then vR = 1.055 * (vR ^ (1 / 2.4)) - 0.055 Else vR = vR * 12.92
            If vG > 0.0031308 Then vG = 1.055 * (vG ^ (1 / 2.4)) - 0.055 Else vG = vG * 12.92
            If vB > 0.0031308 Then vB = 1.055 * (vB ^ (1 / 2.4)) - 0.055 Else vB = vB * 12.92

            Result.Red = Math.Max(Math.Min(vR, 1), 0)
            Result.Green = Math.Max(Math.Min(vG, 1), 0)
            Result.Blue = Math.Max(Math.Min(vB, 1), 0)

            Return Result
        End Function

        Public Shared Function XYZ2HSV(ByVal X As Double, ByVal Y As Double, ByVal Z As Double) As HSVTriplet
            Return RGB2HSV(XYZ2RGB(X, Y, Z))
        End Function

        Public Shared Function XYZ2HSV(ByVal XYZ As XYZTriplet) As HSVTriplet
            Return RGB2HSV(XYZ2RGB(XYZ.X, XYZ.Y, XYZ.Z))
        End Function

        Public Shared Function XYZ2Lab(ByVal X As Double, ByVal Y As Double, ByVal Z As Double) As LabTriplet
            Return XYZ2Lab(X, Y, Z, WhitePointTransformer.StandardWP_D65)
        End Function

        Public Shared Function XYZ2Lab(ByVal XYZ As XYZTriplet, ByVal RefWhite As XYZTriplet) As LabTriplet
            Return XYZ2Lab(XYZ.X, XYZ.Y, XYZ.Z, RefWhite)
        End Function

        Public Shared Function XYZ2Lab(ByVal XYZ As XYZTriplet) As LabTriplet
            Return XYZ2Lab(XYZ, WhitePointTransformer.StandardWP_D65)
        End Function

        Public Shared Function XYZ2Lab(ByVal X As Double, ByVal Y As Double, ByVal Z As Double, ByVal RefWhite As XYZTriplet) As LabTriplet
            Const epsilon As Double = 216.0 / 24389.0
            Const Kappa As Double = 24389.0 / 27.0
            Dim Lab As LabTriplet
            Dim vX As Double = X '/ RefWhite.X
            Dim vY As Double = Y '/ RefWhite.Y
            Dim vZ As Double = Z '/ RefWhite.Z

            If (vX > epsilon) Then vX = vX ^ (1 / 3) Else vX = (Kappa * vX + 16) / 116
            If (vY > epsilon) Then vY = vY ^ (1 / 3) Else vY = (Kappa * vY + 16) / 116
            If (vZ > epsilon) Then vZ = vZ ^ (1 / 3) Else vZ = (Kappa * vZ + 16) / 116
            Lab.Luminance = Math.Min(Math.Max((1.16 * vY) - 0.16, 0), 1)
            Lab.a = Math.Min(Math.Max(5 * (vX - vY), -1.2), 1.2)
            Lab.b = Math.Min(Math.Max(2 * (vY - vZ), -1.2), 1.2)

            Return Lab
        End Function


        Public Shared Function XYZ2Luv(ByVal XYZ As XYZTriplet, ByVal RefWhite As XYZTriplet) As LuvTriplet
            Dim Luv As LuvTriplet
            Dim vU As Double = (4 * XYZ.X) / (XYZ.X + 15 * XYZ.Y + 3 * XYZ.Z)
            Dim vV As Double = (4 * XYZ.Y) / (XYZ.X + 15 * XYZ.Y + 3 * XYZ.Z)

            Dim refU As Double = (4 * RefWhite.X) / (RefWhite.X + 15 * RefWhite.Y + 3 * RefWhite.Z)
            Dim refV As Double = (4 * RefWhite.Y) / (RefWhite.X + 15 * RefWhite.Y + 3 * RefWhite.Z)

            Luv.Luminance = Math.Min(Math.Max(1.16 * Math.Pow(XYZ.Y / RefWhite.Y, 1 / 3) - 0.16, 0), 1)
            Luv.u = Math.Min(Math.Max(13 * Luv.Luminance * (vU - refU), 0), 1)
            Luv.v = Math.Min(Math.Max(13 * Luv.Luminance * (vV - refV), 0), 1)

            Return Luv
        End Function

        Public Shared Function XYZ2xyY(ByVal XYZ As XYZTriplet) As xyYTriplet
            Dim SumXYZ As Double = XYZ.X + XYZ.Y + XYZ.Z
            Return New xyYTriplet(XYZ.X / SumXYZ, XYZ.Y / SumXYZ, XYZ.Y)
        End Function

#End Region

#Region "From HSV"

        Public Shared Function HSV2RGB(ByVal HSV As HSVTriplet) As RGBTriplet
            Dim RGB As RGBTriplet

            If HSV.Saturation > 0 Then
                Dim vH As Double = HSV.Hue * 6

                If (vH >= 6) Then vH = 0
                Dim i As Integer = Math.Floor(vH)
                Dim v1 As Double = HSV.Value * (1 - HSV.Saturation)
                Dim v2 As Double = HSV.Value * (1 - HSV.Saturation * (vH - i))
                Dim v3 As Double = HSV.Value * (1 - HSV.Saturation * (1 - (vH - i)))
                If (i = 0) Then
                    RGB.Red = HSV.Value : RGB.Green = v3 : RGB.Blue = v1
                ElseIf (i = 1) Then
                    RGB.Red = v2 : RGB.Green = HSV.Value : RGB.Blue = v1
                ElseIf (i = 2) Then
                    RGB.Red = v1 : RGB.Green = HSV.Value : RGB.Blue = v3
                ElseIf (i = 3) Then
                    RGB.Red = v1 : RGB.Green = v2 : RGB.Blue = HSV.Value
                ElseIf (i = 4) Then
                    RGB.Red = v3 : RGB.Green = v1 : RGB.Blue = HSV.Value
                Else
                    RGB.Red = HSV.Value : RGB.Green = v1 : RGB.Blue = v2
                End If
            Else
                RGB.Red = HSV.Value
                RGB.Green = HSV.Value
                RGB.Blue = HSV.Value
            End If
            Return RGB
        End Function

        Public Shared Function HSV2XYZ(ByVal HSV As HSVTriplet) As XYZTriplet
            Return RGB2XYZ(HSV2RGB(HSV))
        End Function

        Public Shared Function HSV2Lab(ByVal HSV As HSVTriplet) As LabTriplet
            Return XYZ2Lab(HSV2XYZ(HSV))
        End Function

#End Region

        Public Shared Function xyY2XYZ(ByVal XYy As xyYTriplet) As XYZTriplet
            If XYy.XYZ_Y = 0 Then Return New XYZTriplet(0, 0, 0)

            Dim XYZ As XYZTriplet
            XYZ.X = XYy.x * XYy.XYZ_Y / XYy.y
            XYZ.Y = XYy.XYZ_Y
            XYZ.Z = (1 - XYy.x - XYy.y) * XYy.XYZ_Y / XYy.y
            Return XYZ
        End Function

#Region "From Lab"

        Public Shared Function Lab2RGB(ByVal Lab As LabTriplet) As RGBTriplet
            Return XYZ2RGB(Lab2XYZ(Lab))
        End Function

        Public Shared Function Lab2XYZ(ByVal Lab As LabTriplet) As XYZTriplet
            Return Lab2XYZ(Lab, WhitePointTransformer.StandardWP_D65)
        End Function

        Public Shared Function Lab2XYZ(ByVal Lab As LabTriplet, ByVal RefWhite As XYZTriplet) As XYZTriplet
            Const epsilon As Double = 216.0 / 24389.0
            Const Kappa As Double = 24389.0 / 27.0

            Dim vY As Double = (Lab.Luminance + 0.16) / 1.16
            Dim vX As Double = Lab.a / 5 + vY
            Dim vZ As Double = vY - Lab.b / 2

            If (vX ^ 3 > epsilon) Then vX = vX ^ 3 Else vX = (116 * vX - 16) / Kappa
            If (Lab.Luminance * 100 > epsilon * Kappa) Then vY = vY ^ 3 Else vY = Lab.Luminance / Kappa
            If (vZ ^ 3 > epsilon) Then vZ = vZ ^ 3 Else vZ = (116 * vZ - 16) / Kappa


            'Return New XYZTriplet(vX * RefWhite.X, vY * RefWhite.Y, vZ * RefWhite.Z)
            Return New XYZTriplet(vX, vY, vZ)
        End Function

#End Region


        Public Shared Function RGBE2RGB(ByVal RGBE As RGBEQuadlet) As RGBTriplet
            Return RGBE2RGB(RGBE.Red, RGBE.Green, RGBE.Blue, RGBE.Exponent)
        End Function

        Public Shared Function RGBE2RGB(ByVal R As Byte, ByVal G As Byte, ByVal B As Byte, ByVal E As Byte) As RGBTriplet
            Dim f As Double
            Dim result As Imaging.RGBTriplet
            If E <> 0 Then 'nonzero pixel
                f = ANSI_C_Math.ldexp(1.0, CInt(E) - (128 + 8))
                result.Red = R * f
                result.Green = G * f
                result.Blue = B * f
            End If

            Return result
        End Function


        Public Shared Function LabDelta(ByVal LabColor1 As LabTriplet, ByVal LabColor2 As LabTriplet) As Double
            Return Math.Sqrt((LabColor1.Luminance - LabColor2.Luminance) ^ 2 + _
                             (LabColor1.a - LabColor2.a) ^ 2 + _
                             (LabColor1.b - LabColor2.b) ^ 2)
        End Function

#End Region


    End Class

End Namespace